%%bash
wget http://kodomo.cmm.msu.ru/~golovin/bilayer/dppc.itp
wget http://kodomo.cmm.msu.ru/~golovin/bilayer/lipid.itp
wget http://kodomo.cmm.msu.ru/~golovin/bilayer/dppc.gro
wget http://kodomo.cmm.msu.ru/~golovin/bilayer/b.top 
wget http://kodomo.cmm.msu.ru/~golovin/bilayer/em.mdp
wget http://kodomo.cmm.msu.ru/~golovin/bilayer/pr.mdp
wget http://kodomo.cmm.msu.ru/~golovin/bilayer/md.mdp

! genconf -f dppc.gro -o b_64.gro -nbox 4 4 4

! editconf -f dppc.gro -o dppc.pdb
! editconf -f b_64.gro -o b_64.pdb

Посмотрим в PyMol

Посмотри 64 молекулы

! editconf -f b_64.gro -o b_ec -d 0.5 

! grompp -f em -c b_ec -p b -o b_em -maxwarn 2
! mdrun -deffnm b_em -v

Посмотрим изменение силы

Fmax=4.37970e+05

Fmax=6.16887e+02

! genbox -cp b_em -p b -cs spc216 -o b_s

! grompp -f pr -c b_s -p b -o b_pr -maxwarn 1
! mdrun -deffnm b_pr -v

! grompp -f em -c b_s -p b -o b_empr -maxwarn 1
! mdrun -deffnm b_empr -v

! grompp -f pr -c b_empr -p b -o b_pr -maxwarn 1
!mdrun -deffnm b_pr -v

Переформатируем в pdb

! editconf -f b_pr.gro -o b_pr.pdb
! editconf -f b_s.gro -o b_s.pdb

Посмотрим b_pr

И теперь b_s

В целом, существенной разницы нет, возможно, есть небольшая разница в количество и положениях молекул воды.

Теперь посмотрим результат моделирования

! echo 2 | trjconv -f b_md.xtc -s b_md.tpr -o b_pbc_1.pdb -skip 20 

Попробуем второй способ визуализации

! echo 2 | trjconv -f b_md.xtc -s b_md.tpr -o b_pbc_1.pdb -skip 20 -pbc mol

Посмотрим, что получилось

from IPython.display import Video
Video('b_pbc_1.mpg')

Видим, что довольно быстро, уже в районе 20 кадра образуется бислой.

! echo 0 | g_traj -f b_md.xtc -s b_md.tpr -ob box_1.xvg

Прочитаем результат и построим результат с учётом того, что нормалью к поверхности является ось X

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
%matplotlib inline

data = np.loadtxt('box_1.xvg', comments=['#', '@'])[:,:4]

plt.plot(data[:, 0], data[:, 1] * data[:, 3] / 64)
plt.title('square')
plt.xlabel('time')

Text(0.5,0,u'time')

Как видно по графику, нормализованная площадь падает со временем, это значит, что действительно происходит сборка липидного слоя.

Теперь построим график изменения гидрофобной и гидрофильной поверхности

! g_sas -f b_md.xtc -s b_md.tpr -o sas_b.xvg -dt 100

data = np.loadtxt('sas_b.xvg', comments=['#', '@'])
time = data[:,0]

plt.plot(time, data[:,1], color='green')
plt.plot(time, data[:,2], color='blue')
patch1 = mpatches.Patch(color='green', label='Hydrophobic')
patch2 = mpatches.Patch(color='blue', label='Hydrophilic')
plt.legend(handles=[patch1,patch2])
plt.title('Hydrophobic and hydrophilic surface square')
plt.xlabel('Time')

Text(0.5,0,u'Time')

Видим, что площадь гидрофобных поверхностей падает быстрее, чем площадь гидрофильных, это ожидаемо, так как в процессе сборки липидного слоя молекулы поворачиваются гидрофобными частями внутрь слоя, так как контакты гидробочных поверхностей с молекулами воды энергетически невыгодны. Площадь гидрофильных поверхностей тоже немного падает за счёт компактификации струкруры.

Теперь посмотрим на изменение меры порядка.

!wget http://kodomo.cmm.msu.ru/~golovin/bilayer/sn1.ndx 

! g_order -s b_md -f b_md.xtc -o ord_end.xvg -n sn1.ndx -b 45000 -d X
! g_order -s b_md -f b_md.xtc -o ord_start.xvg -n sn1.ndx -e 5000 -d X

import pandas as pd
start = pd.read_csv('ord_start.xvg', header = None, skiprows = 12, sep='\s+')
end = pd.read_csv('ord_end.xvg', header = None, skiprows = 12, sep='\s+')

plt.plot(start[0], start[3], marker='.', label='start')
plt.plot(end[0], end[3], marker='.', label='end')
plt.xlabel('atom')
plt.ylabel('order')
plt.legend()

<matplotlib.legend.Legend at 0x7fd1906ad910>

Видим, что порядок в результате сборки в конце увеличилась мера порядка, это значит, что в результате сборки структура стала более упорядоченной.